R Notebook

Setup

Run AlphaDiversity in scratchnotebooks

#source(here::here("RScripts", "InitialProcessing_3.R"))
source(here::here("RLibraries", "ChesapeakePersonalLibrary.R"))

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✓ ggplot2 3.3.5     ✓ purrr   0.3.4
✓ tibble  3.1.6     ✓ dplyr   1.0.7
✓ tidyr   1.1.4     ✓ stringr 1.4.0
✓ readr   2.0.0     ✓ forcats 0.5.1
── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

ksource(here::here("ActiveNotebooks", "AlphaDiversity.Rmd"))

processing file: /home/jacob/Projects/ChesapeaakeMainstemAnalysis/ActiveNotebooks/AlphaDiversity.Rmd

  |                                                                                                                                                                                       
  |                                                                                                                                                                                 |   0%
  |                                                                                                                                                                                       
  |......                                                                                                                                                                           |   3%
  |                                                                                                                                                                                       
  |............                                                                                                                                                                     |   7%
  |                                                                                                                                                                                       
  |..................                                                                                                                                                               |  10%
  |                                                                                                                                                                                       
  |........................                                                                                                                                                         |  13%
  |                                                                                                                                                                                       
  |..............................                                                                                                                                                   |  17%
  |                                                                                                                                                                                       
  |...................................                                                                                                                                              |  20%
  |                                                                                                                                                                                       
  |.........................................                                                                                                                                        |  23%
  |                                                                                                                                                                                       
  |...............................................                                                                                                                                  |  27%
  |                                                                                                                                                                                       
  |.....................................................                                                                                                                            |  30%
  |                                                                                                                                                                                       
  |...........................................................                                                                                                                      |  33%
  |                                                                                                                                                                                       
  |.................................................................                                                                                                                |  37%
  |                                                                                                                                                                                       
  |.......................................................................                                                                                                          |  40%
  |                                                                                                                                                                                       
  |.............................................................................                                                                                                    |  43%
  |                                                                                                                                                                                       
  |...................................................................................                                                                                              |  47%
  |                                                                                                                                                                                       
  |........................................................................................                                                                                         |  50%
  |                                                                                                                                                                                       
  |..............................................................................................                                                                                   |  53%
  |                                                                                                                                                                                       
  |....................................................................................................                                                                             |  57%
  |                                                                                                                                                                                       
  |..........................................................................................................                                                                       |  60%
  |                                                                                                                                                                                       
  |................................................................................................................                                                                 |  63%
  |                                                                                                                                                                                       
  |......................................................................................................................                                                           |  67%
  |                                                                                                                                                                                       
  |............................................................................................................................                                                     |  70%
  |                                                                                                                                                                                       
  |..................................................................................................................................                                               |  73%
  |                                                                                                                                                                                       
  |........................................................................................................................................                                         |  77%
  |                                                                                                                                                                                       
  |..............................................................................................................................................                                   |  80%
  |                                                                                                                                                                                       
  |....................................................................................................................................................                             |  83%
  |                                                                                                                                                                                       
  |.........................................................................................................................................................                        |  87%
  |                                                                                                                                                                                       
  |...............................................................................................................................................................                  |  90%
  |                                                                                                                                                                                       
  |.....................................................................................................................................................................            |  93%
  |                                                                                                                                                                                       
  |...........................................................................................................................................................................      |  97%
  |                                                                                                                                                                                       
  |.................................................................................................................................................................................| 100%

output file: /tmp/RtmpWB71Zc/file512e51983b76

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

Attaching package: ‘flextable’

The following object is masked from ‘package:purrr’:

    compose

Warning: Assuming taxa are rows

library(vegan)

Loading required package: permute
Loading required package: lattice
This is vegan 2.6-3

library(cowplot)
library(flextable)
library(ftExtra)

This file is dedicated to conventional, non div-net/breakaway stats, since breakaway seems to choke on this data.

Reshape back into an ASV matrix, but this time correcting for total abundance

nonSpikes %>% head

raDf <- nonSpikes %>% pivot_wider(id_cols = ID, names_from = ASV, values_from = RA, values_fill = 0)
raMat <- raDf %>% column_to_rownames("ID")

raMat1 <- raMat %>% as.matrix()

countMat <-  nonSpikes %>%
  pivot_wider(id_cols = ID, names_from = ASV, values_from = reads, values_fill = 0) %>%
  column_to_rownames("ID") %>% as.matrix()

seqDep <- countMat %>% apply(1, sum)
min(seqDep)

[1] 852

sampleRichness <- rarefy(countMat, min(seqDep))

rarefy everything to the minimum depth (852)

countRare <- rrarefy(countMat, min(seqDep))

Gamma diversity

specpool(countRare)

Doesn’t finish

#specpool(countMat)

Calculate diversity indeces

All richness estimates

richnessRare <- estimateR(countRare)

Shannon diversity

shan <- diversity(countRare)
shan

 3-1-B-0-2  3-1-B-1-2  3-1-B-180   3-1-B-20    3-1-B-5  3-1-B-500   3-1-B-53  3-1-S-0-2  3-1-S-1-2  3-1-S-180   3-1-S-20    3-1-S-5  3-2-B-0-2  3-2-B-1-2  3-2-B-180   3-2-B-20    3-2-B-5 
  4.266589   5.066675   4.748195   5.882121   5.213130   3.866813   5.519614   4.504019   4.997699   4.689925   5.265848   4.805358   4.315456   4.623881   4.559119   5.040415   5.448241 
 3-2-B-500   3-2-B-53  3-2-S-0-2  3-2-S-1-2  3-2-S-180   3-2-S-20    3-2-S-5  3-2-S-500   3-2-S-53  3-3-B-0-2  3-3-B-1-2  3-3-B-180   3-3-B-20    3-3-B-5  3-3-B-500   3-3-B-53  3-3-S-180 
  4.970669   4.754482   3.681555   4.892022   4.667467   5.000035   5.274235   4.862903   4.293725   4.391800   4.959762   3.460213   5.726654   5.328888   5.316755   5.443589   5.054444 
  3-3-S-20  3-3-S-500   3-3-S-53  4-3-B-0-2  4-3-B-1-2  4-3-B-180   4-3-B-20    4-3-B-5  4-3-B-500   4-3-B-53  4-3-O-1-2  4-3-O-180    4-3-O-5  4-3-O-500   4-3-O-53  4-3-S-0-2  4-3-S-180 
  4.906939   4.885977   4.377182   4.224096   4.916432   4.431618   4.408475   4.487108   4.300986   4.176164   4.968473   4.647133   5.112024   3.975860   4.573170   2.805419   4.498084 
  4-3-S-20  4-3-S-500   4-3-S-53  5-1-S-1-2  5-1-S-180   5-1-S-20    5-1-S-5  5-1-S-500   5-1-S-53  5-5-B-0-2  5-5-B-180    5-5-B-5  5-5-B-500   5-5-B-53  5-5-S-180    5-5-S-5  5-5-S-500 
  4.633447   4.763641   4.366296   4.403596   4.477610   4.171291   3.909744   4.041151   4.134044   4.665525   5.098082   5.398357   5.076910   5.036173   4.295120   4.994707   4.830118 
  5-5-S-53 C_5P1B_0P2 C_5P1B_180 C_5P1B_1P2  C_5P1B_20 C_5P1B_500  C_5P1B_53 
  4.247241   4.407919   4.839992   4.888112   5.396460   4.679991   5.232058 

Evenness

pielouJ <- shan/richnessRare["S.chao1",]
pielouJ

  3-1-B-0-2   3-1-B-1-2   3-1-B-180    3-1-B-20     3-1-B-5   3-1-B-500    3-1-B-53   3-1-S-0-2   3-1-S-1-2   3-1-S-180    3-1-S-20     3-1-S-5   3-2-B-0-2   3-2-B-1-2   3-2-B-180 
0.010796521 0.008379501 0.012230611 0.003107806 0.006670746 0.042030574 0.008948671 0.011785878 0.006815941 0.008871841 0.007074017 0.007942428 0.013081311 0.011009240 0.012591089 
   3-2-B-20     3-2-B-5   3-2-B-500    3-2-B-53   3-2-S-0-2   3-2-S-1-2   3-2-S-180    3-2-S-20     3-2-S-5   3-2-S-500    3-2-S-53   3-3-B-0-2   3-3-B-1-2   3-3-B-180    3-3-B-20 
0.006497564 0.006502422 0.007438591 0.008844683 0.019571778 0.006893855 0.008818684 0.008837454 0.009265714 0.009464584 0.014658199 0.014525960 0.006912331 0.062066595 0.005167280 
    3-3-B-5   3-3-B-500    3-3-B-53   3-3-S-180    3-3-S-20   3-3-S-500    3-3-S-53   4-3-B-0-2   4-3-B-1-2   4-3-B-180    4-3-B-20     4-3-B-5   4-3-B-500    4-3-B-53   4-3-O-1-2 
0.006820821 0.007008356 0.006606406 0.010298444 0.010244131 0.010936137 0.011279040 0.012088416 0.009557943 0.011023255 0.009479392 0.008788068 0.010645362 0.009541883 0.008057076 
  4-3-O-180     4-3-O-5   4-3-O-500    4-3-O-53   4-3-S-0-2   4-3-S-180    4-3-S-20   4-3-S-500    4-3-S-53   5-1-S-1-2   5-1-S-180    5-1-S-20     5-1-S-5   5-1-S-500    5-1-S-53 
0.010062029 0.008070617 0.015028669 0.009695583 0.124685269 0.010941017 0.008496069 0.010528326 0.013690750 0.010095361 0.012491681 0.013629647 0.017070925 0.014583301 0.015952761 
  5-5-B-0-2   5-5-B-180     5-5-B-5   5-5-B-500    5-5-B-53   5-5-S-180     5-5-S-5   5-5-S-500    5-5-S-53  C_5P1B_0P2  C_5P1B_180  C_5P1B_1P2   C_5P1B_20  C_5P1B_500   C_5P1B_53 
0.008436353 0.007538475 0.006170381 0.011690997 0.009564602 0.014591206 0.009765933 0.009603814 0.014900732 0.006386277 0.007006523 0.006184548 0.003893466 0.008552268 0.005238606 

Combine diversity data

diversityData <- sampleData %>% 
  left_join(richnessRare %>% t() %>% as.data.frame() %>% rownames_to_column("ID"), by = "ID") %>%
  left_join(shan %>% enframe(name = "ID", value = "shannonH"), by = "ID") %>%
  left_join(pielouJ %>% enframe(name = "ID", value = "pielouJ"), by = "ID") %>%
  arrange(Size_Class)

Generate plots of diversity estimates

Parameters for all plots

plotSpecs <- list(
  facet_wrap(~Depth, ncol = 1) ,
  theme_bw(base_size = 16) ,
  geom_point(size = 4) ,
  geom_path(aes(color = as.factor(Station))) ,
  scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
  #scale_y_log10nice() ,
  scale_shape_manual(values = rep(21:25, 2)) ,
  scale_fill_viridis_d(option = "plasma") ,
  scale_color_viridis_d(option = "plasma") ,
  labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
  theme(legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 90, vjust = .5),
        axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)

Observed species counts, on rarefied data

plotObs <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = S.obs, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +ylab("Observed ASVs (Rarefied)")#+ scale_y_log10()
plotObs

Estemated Chao1 Richness

plotChao1 <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = S.chao1, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  geom_errorbar(aes(ymin = S.chao1 -2 * se.chao1, ymax = S.chao1 + 2* se.chao1), width = -.1) + 
  scale_y_log10() +
  ylab("Richness (Chao1)")
plotChao1

Shannon diversity

plotShan <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = shannonH, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  ylab("Diversity (Shannon H)") +
  lims(y = c(2.5, 6))
plotShan

Evenness

plotPielou <- diversityData %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +scale_y_log10() +ylab("Evenness (PielouJ)")
plotPielou

All plots together

plotAlpha <- plot_grid(plotObs, plotChao1, plotShan, plotPielou, nrow = 1, labels = LETTERS)
plotAlpha

ggsave(here::here("Figures", "ConventionalAlpha.png"), plotAlpha, width = 11, height = 4)

Do we see trends with lat and size?

Observed Species

Rarefied

obsMod <- lm(S.obs ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(obsMod)

Call:
lm(formula = S.obs ~ log(Size_Class) + I(log(Size_Class)^2) + 
    I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)

Residuals:
     Min       1Q   Median       3Q      Max 
-222.130  -40.254    0.968   47.519  200.464 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)          294642.461 134921.667   2.184 0.032379 *  
log(Size_Class)          33.747      8.174   4.128 0.000101 ***
I(log(Size_Class)^2)     -6.605      1.521  -4.343 4.72e-05 ***
lat                  -15345.129   7030.594  -2.183 0.032470 *  
I(lat^2)                199.852     91.525   2.184 0.032397 *  
depth                     5.038      3.238   1.556 0.124293    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 76.85 on 69 degrees of freedom
Multiple R-squared:  0.2783,    Adjusted R-squared:  0.226 
F-statistic: 5.322 on 5 and 69 DF,  p-value: 0.0003409

Richness

chao1Mod <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(chao1Mod)

Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) + 
    I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)

Residuals:
    Min      1Q  Median      3Q     Max 
-516.16 -139.22  -14.04  140.31 1146.14 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)          952088.842 448528.950   2.123 0.037370 *  
log(Size_Class)          85.540     27.175   3.148 0.002430 ** 
I(log(Size_Class)^2)    -18.549      5.056  -3.669 0.000476 ***
lat                  -49615.130  23372.263  -2.123 0.037359 *  
I(lat^2)                646.322    304.262   2.124 0.037237 *  
depth                    19.787     10.763   1.838 0.070303 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 255.5 on 69 degrees of freedom
Multiple R-squared:  0.2193,    Adjusted R-squared:  0.1627 
F-statistic: 3.877 on 5 and 69 DF,  p-value: 0.003731

chao1ModSimple <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2), data = diversityData)
summary(chao1ModSimple)

Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2), 
    data = diversityData)

Residuals:
    Min      1Q  Median      3Q     Max 
-451.11 -148.26  -32.91  126.61 1234.57 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)           569.961     46.362  12.294  < 2e-16 ***
log(Size_Class)        86.118     27.533   3.128 0.002542 ** 
I(log(Size_Class)^2)  -18.923      5.122  -3.695 0.000426 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 259.2 on 72 degrees of freedom
Multiple R-squared:  0.1617,    Adjusted R-squared:  0.1384 
F-statistic: 6.943 on 2 and 72 DF,  p-value: 0.001748

Shannon Diversity

shanMod <- lm(shannonH ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)

summary(shanMod)

Call:
lm(formula = shannonH ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityData)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.32009 -0.22746  0.04123  0.30999  0.76552 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           1.887e+03  7.925e+02   2.382   0.0200 *  
log(Size_Class)       2.051e-01  4.801e-02   4.273 6.06e-05 ***
I(log(Size_Class)^2) -3.764e-02  8.933e-03  -4.213 7.48e-05 ***
lat                  -9.807e+01  4.129e+01  -2.375   0.0203 *  
I(lat^2)              1.277e+00  5.376e-01   2.375   0.0204 *  
depth                 2.466e-02  1.902e-02   1.297   0.1991    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4514 on 69 degrees of freedom
Multiple R-squared:  0.3087,    Adjusted R-squared:  0.2586 
F-statistic: 6.161 on 5 and 69 DF,  p-value: 8.916e-05

Evenness

pielouMod <- lm(pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(pielouMod)

Call:
lm(formula = pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityData)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.016223 -0.005097 -0.002264  0.000799  0.099586 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)          -2.365e+01  2.621e+01  -0.902   0.3700  
log(Size_Class)      -4.048e-03  1.588e-03  -2.549   0.0130 *
I(log(Size_Class)^2)  7.103e-04  2.955e-04   2.404   0.0189 *
lat                   1.232e+00  1.366e+00   0.902   0.3704  
I(lat^2)             -1.601e-02  1.778e-02  -0.901   0.3709  
depth                -3.282e-04  6.290e-04  -0.522   0.6034  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01493 on 69 degrees of freedom
Multiple R-squared:  0.1032,    Adjusted R-squared:  0.03826 
F-statistic: 1.589 on 5 and 69 DF,  p-value: 0.1748

uomisto H (2010a). “A diversity of beta diver- sities: straightening up a concept gone awry. 1. Defining beta diversity as a function of alpha and gamma diversity.” Ecography, 33, 2–2

Prediction plots

Observed Species

predict(obsMod, se.fit = TRUE)

$fit
       1        2        3        4        5        6        7        8        9       10       11       12       13       14       15       16       17       18       19       20 
225.2948 225.2948 191.7909 191.7909 205.5804 152.9884 152.9884 184.9565 198.0723 302.6501 302.6501 269.1462 269.1462 282.9357 230.3437 230.3437 262.3118 262.3118 333.9210 333.9210 
      21       22       23       24       25       26       27       28       29       30       31       32       33       34       35       36       37       38       39       40 
300.4171 300.4171 314.2066 261.6146 261.6146 293.5827 306.6985 306.6985 338.5363 338.5363 305.0324 305.0324 318.8220 318.8220 266.2299 266.2299 298.1980 298.1980 326.5841 293.0801 
      41       42       43       44       45       46       47       48       49       50       51       52       53       54       55       56       57       58       59       60 
293.0801 306.8697 306.8697 254.2777 254.2777 254.2777 286.2457 286.2457 299.3615 299.3615 293.8447 293.8447 260.3408 260.3408 274.1303 274.1303 221.5383 221.5383 221.5383 253.5064 
      61       62       63       64       65       66       67       68       69       70       71       72       73       74       75 
253.5064 266.6221 266.6221 251.3431 217.8392 217.8392 231.6287 231.6287 179.0367 179.0367 179.0367 211.0048 211.0048 224.1206 224.1206 

$se.fit
 [1] 26.88498 26.88498 26.15198 26.15198 26.01242 27.67418 27.67418 29.22809 33.90698 18.85945 18.85945 18.24545 18.24545 17.13460 19.96313 19.96313 21.42036 21.42036 19.17995 19.17995
[21] 18.71991 18.71991 17.16981 20.00965 20.00965 21.12357 28.06863 28.06863 19.17779 19.17779 18.68694 18.68694 16.89659 16.89659 19.53551 19.53551 20.49627 20.49627 18.49063 17.85404
[41] 17.85404 15.93670 15.93670 18.37221 18.37221 18.37221 19.35941 19.35941 26.45249 26.45249 19.23795 19.23795 18.35438 18.35438 16.61635 16.61635 18.35588 18.35588 18.35588 19.42396
[61] 19.42396 26.01457 26.01457 24.36577 23.41633 23.41633 22.25921 22.25921 23.05203 23.05203 23.05203 24.04762 24.04762 29.11453 29.11453

$df
[1] 69

$residual.scale
[1] 76.8497

diversityData$pred_obs = predict(obsMod, se.fit = TRUE)$fit
diversityData$se_obs = predict(obsMod, se.fit = TRUE)$se.fit

plotSpecs2 <- list(
  facet_wrap(~Depth, ncol = 1) ,
  theme_bw(base_size = 16) ,
  #geom_point(size = 4) ,
  geom_path(aes(color = as.factor(Station))) ,
  scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
  #scale_y_log10nice() ,
  scale_shape_manual(values = rep(21:25, 2)) ,
  scale_fill_viridis_d(option = "plasma") ,
  scale_color_viridis_d(option = "plasma") ,
  labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
  theme(legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 90, vjust = .5),
        axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)

plotObs_pred <-  diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_obs, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_obs - 2 * se_obs, yend = pred_obs + 2 * se_obs, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5)  +
  plotSpecs2 + ylab("Predicted  ASVs") 
plotObs_pred

Richness

predict(chao1Mod, se.fit = TRUE)

$fit
       1        2        3        4        5        6        7        8        9       10       11       12       13       14       15       16       17       18       19       20 
471.0415 471.0415 357.8857 357.8857 442.2027 277.6790 277.6790 403.5903 380.9771 671.7413 671.7413 558.5855 558.5855 642.9025 478.3788 478.3788 604.2901 604.2901 746.3854 746.3854 
      21       22       23       24       25       26       27       28       29       30       31       32       33       34       35       36       37       38       39       40 
633.2296 633.2296 717.5466 553.0229 553.0229 678.9342 656.3210 656.3210 746.5474 746.5474 633.3915 633.3915 717.7085 717.7085 553.1848 553.1848 679.0961 679.0961 703.9826 590.8268 
      41       42       43       44       45       46       47       48       49       50       51       52       53       54       55       56       57       58       59       60 
590.8268 675.1438 675.1438 510.6201 510.6201 510.6201 636.5314 636.5314 613.9182 613.9182 600.7490 600.7490 487.5932 487.5932 571.9102 571.9102 407.3865 407.3865 407.3865 533.2978 
      61       62       63       64       65       66       67       68       69       70       71       72       73       74       75 
533.2978 510.6846 510.6846 471.9554 358.7995 358.7995 443.1165 443.1165 278.5929 278.5929 278.5929 404.5041 404.5041 381.8910 381.8910 

$se.fit
 [1]  89.37550  89.37550  86.93875  86.93875  86.47479  91.99908  91.99908  97.16485 112.71919  62.69571  62.69571  60.65454  60.65454  56.96169  66.36473  66.36473  71.20912  71.20912
[19]  63.76118  63.76118  62.23181  62.23181  57.07873  66.51939  66.51939  70.22248  93.31038  93.31038  63.75398  63.75398  62.12222  62.12222  56.17044  56.17044  64.94318  64.94318
[37]  68.13710  68.13710  61.46961  59.35336  59.35336  52.97941  52.97941  61.07595  61.07595  61.07595  64.35777  64.35777  87.93774  87.93774  63.95397  63.95397  61.01666  61.01666
[55]  55.23882  55.23882  61.02166  61.02166  61.02166  64.57236  64.57236  86.48194  86.48194  81.00074  77.84443  77.84443  73.99777  73.99777  76.63338  76.63338  76.63338  79.94309
[73]  79.94309  96.78733  96.78733

$df
[1] 69

$residual.scale
[1] 255.4765

diversityData$pred_chao1 = predict(chao1Mod, se.fit = TRUE)$fit
diversityData$se_chao1 = predict(chao1Mod, se.fit = TRUE)$se.fit

plotChao1_pred <-  diversityData %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_chao1, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_chao1 - 2 * se_chao1, yend = pred_chao1 + 2 * se_chao1, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5)  +
  plotSpecs2 + ylab("Predictd Richness (Chao1)") + scale_y_log10()
plotChao1_pred

Shannon Diversity

predict(shanMod, se.fit = TRUE)

$fit
       1        2        3        4        5        6        7        8        9       10       11       12       13       14       15       16       17       18       19       20 
4.495853 4.495853 4.294730 4.294730 4.302302 3.967491 3.967491 4.130228 4.372581 4.959672 4.959672 4.758549 4.758549 4.766122 4.431310 4.431310 4.594047 4.594047 5.156199 5.156199 
      21       22       23       24       25       26       27       28       29       30       31       32       33       34       35       36       37       38       39       40 
4.955075 4.955075 4.962648 4.627837 4.627837 4.790574 5.032926 5.032926 5.200306 5.200306 4.999182 4.999182 5.006755 5.006755 4.671944 4.671944 4.834681 4.834681 5.144714 4.943591 
      41       42       43       44       45       46       47       48       49       50       51       52       53       54       55       56       57       58       59       60 
4.943591 4.951164 4.951164 4.616353 4.616353 4.616353 4.779090 4.779090 5.021442 5.021442 4.973856 4.973856 4.772732 4.772732 4.780305 4.780305 4.445494 4.445494 4.445494 4.608231 
      61       62       63       64       65       66       67       68       69       70       71       72       73       74       75 
4.608231 4.850584 4.850584 4.744786 4.543663 4.543663 4.551236 4.551236 4.216425 4.216425 4.216425 4.379161 4.379161 4.621514 4.621514 

$se.fit
 [1] 0.15791178 0.15791178 0.15360644 0.15360644 0.15278670 0.16254721 0.16254721 0.17167427 0.19915623 0.11077297 0.11077297 0.10716658 0.10716658 0.10064191 0.11725554 0.11725554
[17] 0.12581478 0.12581478 0.11265549 0.11265549 0.10995335 0.10995335 0.10084870 0.11752879 0.11752879 0.12407155 0.16486406 0.16486406 0.11264276 0.11264276 0.10975972 0.10975972
[33] 0.09924390 0.09924390 0.11474388 0.11474388 0.12038702 0.12038702 0.10860666 0.10486760 0.10486760 0.09360588 0.09360588 0.10791114 0.10791114 0.10791114 0.11370956 0.11370956
[49] 0.15537149 0.15537149 0.11299611 0.11299611 0.10780639 0.10780639 0.09759788 0.09759788 0.10781521 0.10781521 0.10781521 0.11408871 0.11408871 0.15279933 0.15279933 0.14311495
[65] 0.13753827 0.13753827 0.13074185 0.13074185 0.13539854 0.13539854 0.13539854 0.14124626 0.14124626 0.17100726 0.17100726

$df
[1] 69

$residual.scale
[1] 0.4513849

diversityData$pred_shanH = predict(shanMod, se.fit = TRUE)$fit
diversityData$se_shanH = predict(shanMod, se.fit = TRUE)$se.fit

plotShannonH_pred <- diversityData %>%

 filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_shanH, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_shanH - 2 * se_shanH, yend = pred_shanH + 2 * se_shanH, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"),  alpha = 0.5)  +
  plotSpecs2 + ylab("Predicted Diversity (Shannon H)") #+ scale_y_log10()
plotShannonH_pred

Evenness

predict(pielouMod, se.fit = TRUE)

$fit
          1           2           3           4           5           6           7           8           9          10          11          12          13          14          15 
0.019650593 0.019650593 0.022081017 0.022081017 0.021538847 0.025099422 0.025099422 0.022609276 0.019055594 0.010581192 0.010581192 0.013011616 0.013011616 0.012469446 0.016030021 
         16          17          18          19          20          21          22          23          24          25          26          27          28          29          30 
0.016030021 0.013539875 0.013539875 0.006620266 0.006620266 0.009050690 0.009050690 0.008508520 0.012069095 0.012069095 0.009578950 0.006025267 0.006025267 0.005542827 0.005542827 
         31          32          33          34          35          36          37          38          39          40          41          42          43          44          45 
0.007973250 0.007973250 0.007431081 0.007431081 0.010991656 0.010991656 0.008501510 0.008501510 0.006419563 0.008849987 0.008849987 0.008307817 0.008307817 0.011868392 0.011868392 
         46          47          48          49          50          51          52          53          54          55          56          57          58          59          60 
0.011868392 0.009378246 0.009378246 0.005824564 0.005824564 0.009427605 0.009427605 0.011858029 0.011858029 0.011315859 0.011315859 0.014876434 0.014876434 0.014876434 0.012386288 
         61          62          63          64          65          66          67          68          69          70          71          72          73          74          75 
0.012386288 0.008832606 0.008832606 0.013569678 0.016000102 0.016000102 0.015457932 0.015457932 0.019018507 0.019018507 0.019018507 0.016528362 0.016528362 0.012974679 0.012974679 

$se.fit
 [1] 0.005223138 0.005223138 0.005080733 0.005080733 0.005053619 0.005376461 0.005376461 0.005678351 0.006587352 0.003663961 0.003663961 0.003544674 0.003544674 0.003328862 0.003878380
[16] 0.003878380 0.004161488 0.004161488 0.003726227 0.003726227 0.003636850 0.003636850 0.003335702 0.003887418 0.003887418 0.004103828 0.005453094 0.005453094 0.003725806 0.003725806
[31] 0.003630446 0.003630446 0.003282621 0.003282621 0.003795304 0.003795304 0.003981958 0.003981958 0.003592307 0.003468633 0.003468633 0.003096136 0.003096136 0.003569302 0.003569302
[46] 0.003569302 0.003761092 0.003761092 0.005139115 0.005139115 0.003737494 0.003737494 0.003565837 0.003565837 0.003228177 0.003228177 0.003566129 0.003566129 0.003566129 0.003773633
[61] 0.003773633 0.005054037 0.005054037 0.004733714 0.004549258 0.004549258 0.004324457 0.004324457 0.004478483 0.004478483 0.004478483 0.004671904 0.004671904 0.005656288 0.005656288

$df
[1] 69

$residual.scale
[1] 0.01493014

diversityData$pred_pielouJ = predict(pielouMod, se.fit = TRUE)$fit
diversityData$se_pielouJ = predict(pielouMod, se.fit = TRUE)$se.fit

plot_pielouJ_pred <- diversityData %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_pielouJ - 2 * se_pielouJ, yend = pred_pielouJ + 2 * se_pielouJ, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"))  +
  plotSpecs2 + ylab("Predicted Evenness (Pielou J)") + scale_y_log10()
plot_pielouJ_pred

Combined prediction plot

plotPredictions <- plot_grid(plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred, nrow = 1, labels = LETTERS)

Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (geom_segment).

plotPredictions

ggsave(here::here("Figures", "ConventionalAlphaPredictions.png"), plotPredictions, width = 11, height = 4)

Even combindeder

plot_grid(plotObs, plotChao1, plotShan, plotPielou,
          plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred,
          nrow = 2, labels = LETTERS)

Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (geom_segment).

Combined summary table

alphaSummary <- tibble(
  metric = c("Observed ASVs", "Richness (Chao1)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
  model = list(obsMod, chao1Mod, shanMod, pielouMod)
)

alphaSummary <- alphaSummary %>%
  mutate(df = map(model, ~broom::tidy(summary(.))))

alphaSummary <- alphaSummary %>%
  select(-model) %>%
  unnest(df)

alphaSummary <- alphaSummary %>%
  rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
  mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
         Term = str_replace(Term, "\\^2", "\\^2\\^"),
         Term = str_replace(Term, "depth", "Depth"),
         Term = str_replace(Term, "lat", "Latitude"),
         Term = str_replace(Term, "_", " ")# BOOKMARK!!
         ) %>%
  mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
           reformat_sci()
         ) %>%
  mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
           reformat_sci()
  ) %>%
  mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
  mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
  rename(`Standard\nError` = `Standard Error`) %>%
  identity()

alphaSummary %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>% bold(i = ~ p< 0.05, j = "p") %>% colformat_md() %>% set_table_properties(layout = "autofit") %>%
  align(j = -c(1:2), align = "right")

Now considering breakaway values

richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.))))

Warning: `rename_()` was deprecated in dplyr 0.7.0.
Please use `rename()` instead.

diversityDataWB <- full_join(diversityData,
                             richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.)))),
                             by = c("ID" = "break_sample_names"), suffix = c("", "_break")) %>%
  mutate(pielouJ2 = shannonH/break_estimate) %>%
  identity()

diversityDataWB

pielouMod2 <- lm(pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityDataWB)
summary(pielouMod2)

Call:
lm(formula = pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityDataWB)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.013948 -0.005120 -0.002519  0.000893  0.105997 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)          -1.940e+01  2.756e+01  -0.704   0.4840  
log(Size_Class)      -3.284e-03  1.670e-03  -1.966   0.0533 .
I(log(Size_Class)^2)  5.739e-04  3.107e-04   1.847   0.0690 .
lat                   1.009e+00  1.436e+00   0.702   0.4849  
I(lat^2)             -1.310e-02  1.870e-02  -0.701   0.4858  
depth                -2.331e-04  6.614e-04  -0.352   0.7256  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.0157 on 69 degrees of freedom
Multiple R-squared:  0.06678,   Adjusted R-squared:  -0.0008485 
F-statistic: 0.9875 on 5 and 69 DF,  p-value: 0.4318

Ok. So the narrative makes sense. Alpha diveristy is driven by variability in richness rather than evenness. Why would we see an effect in chao1 but not breakaway? Because chao1 is more driven by abundant stuff that makes the rarification threshold. My first hunch is that chao1 responds to evenness, but actually that shouldn’t have any effect since there is no evenness variability? Or maybe just that breakaway is more variable (because it detects fine level differences in rare species) and that doesn’t map as nicely with overall patterns.

plotBreak <- diversityDataWB %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  #scale_y_log10()+
  ylab("Richness (Breakaway)")
plotBreak

plotPielou2 <- diversityDataWB %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  #scale_y_log10()+
  ylab("Evenness (PielouJ)")
plotPielou2

Redo predictions for good measure

predict(pielouMod2, se.fit = TRUE)

$fit
            1             2             3             4             5             6             7             8             9            10            11            12            13 
 0.0116712168  0.0116712168  0.0135414829  0.0135414829  0.0133652561  0.0159774714  0.0159774714  0.0139770206  0.0099002360  0.0043201718  0.0043201718  0.0061904378  0.0061904378 
           14            15            16            17            18            19            20            21            22            23            24            25            26 
 0.0060142111  0.0086264264  0.0086264264  0.0066259755  0.0066259755  0.0011013110  0.0011013110  0.0029715770  0.0029715770  0.0027953503  0.0054075656  0.0054075656  0.0034071147 
           27            28            29            30            31            32            33            34            35            36            37            38            39 
-0.0006696698 -0.0006696698  0.0002127168  0.0002127168  0.0020829828  0.0020829828  0.0019067561  0.0019067561  0.0045189714  0.0045189714  0.0025185205  0.0025185205  0.0009083804 
           40            41            42            43            44            45            46            47            48            49            50            51            52 
 0.0027786464  0.0027786464  0.0026024197  0.0026024197  0.0052146350  0.0052146350  0.0052146350  0.0032141841  0.0032141841 -0.0008626004 -0.0008626004  0.0033228151  0.0033228151 
           53            54            55            56            57            58            59            60            61            62            63            64            65 
 0.0051930811  0.0051930811  0.0050168543  0.0050168543  0.0076290696  0.0076290696  0.0076290696  0.0056286188  0.0056286188  0.0015518342  0.0015518342  0.0066561189  0.0085263849 
           66            67            68            69            70            71            72            73            74            75 
 0.0085263849  0.0083501582  0.0083501582  0.0109623735  0.0109623735  0.0109623735  0.0089619227  0.0089619227  0.0048851381  0.0048851381 

$se.fit
 [1] 0.005492160 0.005492160 0.005342421 0.005342421 0.005313910 0.005653380 0.005653380 0.005970819 0.006926639 0.003852676 0.003852676 0.003727246 0.003727246 0.003500318 0.004078139
[16] 0.004078139 0.004375829 0.004375829 0.003918150 0.003918150 0.003824170 0.003824170 0.003507510 0.004087643 0.004087643 0.004315200 0.005733960 0.005733960 0.003917707 0.003917707
[31] 0.003817435 0.003817435 0.003451696 0.003451696 0.003990784 0.003990784 0.004187052 0.004187052 0.003777332 0.003647288 0.003647288 0.003255606 0.003255606 0.003753142 0.003753142
[46] 0.003753142 0.003954810 0.003954810 0.005403809 0.005403809 0.003929997 0.003929997 0.003749498 0.003749498 0.003394447 0.003394447 0.003749805 0.003749805 0.003749805 0.003967997
[61] 0.003967997 0.005314350 0.005314350 0.004977528 0.004783571 0.004783571 0.004547192 0.004547192 0.004709152 0.004709152 0.004709152 0.004912535 0.004912535 0.005947620 0.005947620

$df
[1] 69

$residual.scale
[1] 0.01569913

diversityDataWB$pred_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$fit
diversityDataWB$se_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$se.fit

plot_pielouJ2_pred <- diversityDataWB %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_pielouJ2 - 2 * se_pielouJ2, yend = pred_pielouJ2 + 2 * se_pielouJ2, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"))  +
  plotSpecs2 + ylab("Predicted Evenness (Pielou J2)") #+ scale_y_log10()
plot_pielouJ2_pred

Breakaway richness subplots

plotBreakaway <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  geom_errorbar(aes(ymin = break_lower, ymax = break_upper), width = -.1) + 
  scale_y_log10() +
  ylab("Richness (Breakaway)")
plotBreakaway

#predict(breakLm, se.fit = TRUE)
# doesn't work because built with a different data frame

Why are these not smooth curves?!! What if I redo the model, this time with the same data frame

breakLm2 <- lm(break_estimate ~ log(Size_Class) + I(log(Size_Class) ^2) + lat +  I(lat^2) + depth ,data = diversityDataWB)
breakLm2 %>% summary()

Call:
lm(formula = break_estimate ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityDataWB)

Residuals:
    Min      1Q  Median      3Q     Max 
-2974.5 -1191.2  -151.6   599.9  6768.1 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)          7124615.61 3339862.88   2.133   0.0365 *
log(Size_Class)          244.45     202.35   1.208   0.2312  
I(log(Size_Class)^2)     -75.16      37.65  -1.996   0.0498 *
lat                  -370568.38  174035.93  -2.129   0.0368 *
I(lat^2)                4817.28    2265.61   2.126   0.0371 *
depth                    151.10      80.15   1.885   0.0636 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1902 on 69 degrees of freedom
Multiple R-squared:  0.1414,    Adjusted R-squared:  0.0792 
F-statistic: 2.273 on 5 and 69 DF,  p-value: 0.0567

Note the non statistical significance overall

#predict(breakLm2, se.fit = TRUE)
diversityDataWB$pred_break = predict(breakLm2, se.fit = TRUE)$fit
diversityDataWB$se_break = predict(breakLm2, se.fit = TRUE)$se.fit

plot_break_pred <- diversityDataWB %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
#  filter(Station == 4.3) %>%
  ggplot(aes(x = Size_Class, y = pred_break, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_break - 2 * se_break, yend = pred_break + 2 * se_break, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"))  +
  plotSpecs2 + ylab("Predicted Richness (Breakaway -- LM)") #+ scale_y_log10()
plot_break_pred

Rebuilding combined products

plotAlphaWB <- plot_grid(plotBreakaway, plotShan, plotPielou2, nrow = 1, labels = LETTERS)
plotAlphaWB

ggsave(here::here("Figures", "BreakawayAlpha.png"), plotAlpha, width = 11, height = 4)

Summary table I want both breakaway metrics here

bettaTable <- myBet$table %>% 
  as.data.frame() %>%
  rename(estimate = Estimates,
         `std.error` = `Standard Errors`,
         `p.value`=`p-values`
         ) %>%
  mutate(`statistic` = NA) %>%
  rownames_to_column(var = "term") %>%
  select(term, estimate, std.error, statistic, p.value) %>%
  as_tibble()
bettaTable

alphaSummary2 <- tibble(
  metric = c("Richness (Breakaway -- LM)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
  model = list(breakLm, shanMod, pielouMod2)
)
  
alphaSummary2 <- alphaSummary2 %>%
  mutate(df = map(model, ~broom::tidy(summary(.))))

## Add in willis variables

breakawaySummary <- tibble(
  metric = "Richness (Breakaway -- Betta)",
  model = NULL,
  df = list(bettaTable)
)

alphaSummary2 = bind_rows(breakawaySummary, alphaSummary2)

alphaSummary2 <- alphaSummary2 %>%
  select(-model) %>%
  unnest(df)

alphaSummary2 <- alphaSummary2 %>%
  rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
  mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
         Term = str_replace(Term, "\\^2", "\\^2\\^"),
         Term = str_replace(Term, "depth", "Depth"),
         Term = str_replace(Term, "lat", "Latitude"),
         Term = str_replace(Term, "_", " ")# BOOKMARK!!
         ) %>%
  mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
           reformat_sci()
         ) %>%
  mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
           reformat_sci()
  ) %>%
  mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
  mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
  rename(`Standard\nError` = `Standard Error`) %>%
  identity()



alphaSummary2

alphaTable2 <- alphaSummary2 %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>% bold(i = ~ p< 0.05, j = "p") %>% colformat_md() %>% set_table_properties(layout = "autofit") %>%
  align(j = -c(1:2), align = "right")
alphaTable2

alphaTable2 %>% save_as_docx(path = here::here("Tables", "alphaTable2.docx"))

myBet$table

And finally predictions from richness, diversity evenness again.

plot_grid(plot_break_pred,plotShannonH_pred,plot_pielouJ2_pred, nrow = 1, labels = LETTERS)